arXiv:1507.04052vl [physics.comp-ph] 14Jul2015 


Optimizing the geometrical accuracy of curvilinear 

meshes 

Thomas Touiorge f,J \ Jonathan Lambrechts a,b and Jean-Frangois Remacle a 
a Universite catholique de Louvain, Institute of Mechanics, Materials and Civil 
Engineering (iMMC), Batiment Euler, Avenue Georges Lemaitre 4, 

1348 Louvain-la-Neuve, Belgium 

6 Fonds National de la Recherche Scientifique, rue d’Egmond 5, 1000 Bruxelles, Belgium 

July 16, 2015 


Abstract 

This paper presents a method to generate valid high order meshes 
with optimized geometrical accuracy. The high order meshing proce¬ 
dure starts with a linear mesh, that is subsequently curved without 
taking care of the validity of the high order elements. An optimization 
procedure is then used to both untangle invalid elements and opti¬ 
mize the geometrical accuracy of the mesh. Standard measures of the 
distance between curves are considered to evaluate the geometrical 
accuracy in planar two-dimensional meshes, but they prove compu¬ 
tationally too costly for optimization purposes. A fast estimate of 
the geometrical accuracy, based on Taylor expansions of the curves, is 
introduced. An unconstrained optimization procedure based on this 
estimate is shown to yield significant improvements in the geometrical 
accuracy of high order meshes, as measured by the standard Haudorff 
distance between the geometrical model and the mesh. Several exam¬ 
ples illustrate the beneficial impact of this method on CFD solutions, 
with a particular role of the enhanced mesh boundary smoothness. 


1 Introduction 

The development of high-order numerical technologies for engineering anal¬ 
ysis has been underway for many years now. For example, Discontinuous 
Galerkin methods (DGM) have been thoroughly studied in the literature, 
initially in a theoretical context [7], and now from the application point 
of view m ns]. Compared to standard second-order-accurate numerical 
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schemes, high-order methods exhibit superior efficiency in problems with 
high resolution requirements, because they reach the required accuracy with 
much coarser grids. 

However, many contributions have pointed out that the accuracy of these 
methods can be severely hampered by a too crude discretization of the geom¬ 
etry mum- It is now widely accepted that linear geometrical discretiza¬ 
tions may annihilate the benefits of high-order schemes in cases featuring 
curved geometries, that is, in most cases of engineering and scientific inter¬ 
est. 

This problem has motivated the development of methods for the genera¬ 
tion of high-order meshes, in which curvilinear elements are meant to provide 
sufficient geometrical accuracy on the boundary. Elements are then most 
often defined in a Lagrangian manner by a set of high-order nodes. Until 
now, efforts have mostly been targeted at ensuring the validity of the mesh. 
Indeed, the naive approach consisting in simply curving the boundaries of 
a linear mesh to match the geometry often results in tangled elements [2i| . 
The curvature of the boundary must somehow be “propagated” into the do¬ 
main for all elements to be valid. In the case of locally structured meshes, 
such a situation can be avoided by means of an efficient isoparametric tech¬ 
nique m- For unstructured meshes, untangling procedures based on topo¬ 
logical operations [Si, [T6- '2(l!, mechanical analogies ESHHH or optimization 
procedures mm have been proposed. 

Although the improved representation of the geometry of the domain 
is the prime motivation for the use of high-order meshes, only few authors 
have taken into consideration the quality of the geometrical approximation. 
In the literature, the limited work dedicated to this topic has focused on 
placing adequately the high-order nodes when curving the mesh bound¬ 
aries. Simple techniques include interpolating them between the first-order 
boundary nodes in the parametric space describing the corresponding CAD 
entity [ ; 8J [lOj or projecting them on the geometry from their location on 
the straight-sided element. More sophisticated procedures have also been 
proposed. In Ref. [26], the high-order nodes on boundary edges are interpo¬ 
lated in the physical space through a numerical procedure involving either 
the CAD parametrization (in the case of a mesh edge assigned to an edge 
of the geometric model), or an approximation of the geodesic connecting 
the two first-order vertices (in the case of an edge located on a 3D sur¬ 
face). Nodes located within surface elements are obtained through a more 
sophisticated version of this procedure. Instead of interpolating, Sherwin 
and Peiro [2T] use a mechanical analogy with chains of springs in equilib¬ 
rium that yields the adequate node distribution along geometric curves and 
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geodesics for edge nodes. Two-dimensional nets of springs provide the ap¬ 
propriate distribution of surface element nodes. 

This paper presents a method that makes it possible to build geomet¬ 
rically accurate curvilinear meshes. Unlike previous work reported in the 
literature, the representation of the model by the mesh is formally assessed 
by measuring distances between the geometric model and the corresponding 
high-order mesh boundary, or by evaluating a fast estimate of the geomet¬ 
rical error. The aim of the method is to minimize this geometrical error 
through the use of standard optimization algorithms. Although most of the 
paper deals with two-dimensional meshes, it is shown that the approach can 
easily be extended to three spatial dimensions. 

Consider a model entity C and the mesh entity C m that is meant to 
approximate C. The first questions that arise are how to define a proper 
distance d(C,C m ) between C and C m , and how to compute this distance 
efficiently. Two main definitions for such a distance have been proposed in 
the computational geometry literature, namely the Frechet distance and the 
Hausdorff distance. 

In this context of curvilinear meshing, distances d(C,C m ) that are com¬ 
puted are usually small in comparison with the typical dimension of either 
C or Cm- Consequently, d(C,C m ) has to be computed with high accuracy. In 
this paper, we show that computing standard distances between the mesh 
and the geometry may be too expensive for practical computations. Al¬ 
ternative measures of distance are presented, that are both fast enough to 
compute and sufficiently accurate. An optimization procedure is then devel¬ 
oped to drastically reduce the model-to-mesh distance while enforcing the 
mesh validity. 

The paper is organized as follows. In Section [2j the problem of defining 
and computing a proper model-to-mesh distance is examined. The mesh 
optimization procedure is described in Section [3| Section [f] illustrates the 
method with examples, and the extension of the approach to three dimen¬ 
sions is presented in Section [5j Conclusions are drawn in Section [6j 

2 Model-to-Mesh Distance 

2.1 Distance Between Curves 

2.1.1 Setup 

Consider the following planar parametric curve 

C = {rje [r/ 0 , rjp\ *(r/) € R 2 } 
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and the following p + 1 successive points on C 


Xi = x(rn), with r/o < m < m ■ ■ ■ < V P -1 < dp- 

A curvilinear mesh edge C m is defined as the Lagrange approximation of C 
at order p 

Cm = /£ £ [0,1] H > x m (£) = ^2 £-i P \0 S-R 2 !- (1) 

In ([I]), c[ p \^) is the ith Lagrange polynomial of order p. 

Curves C and C rn that are both bounded by the vertices xq and x p and 
coincide at least at the p + 1 Lagrange points Xi, i = 0,... ,p (see Figure 

El- 



Figure 1: Typical setup: a model edge C and a quadratic mesh edge C m . 


2.1.2 Formal Definitions of Distance 

Define a(t) (resp. /3(t)) as an arbitrary continuous nondecreasing function 
from t € [0,1] onto rj E [770 ,Vp\ (resp. ^ E [0,1]). The Frechet distance 
between C and C m is defined as 

d F (C,Cm ) = inf max || x m (/3(t)) -x(a(t))||. 
q,/ 3 te[o,i] 

There is a standard interpretation of the Frechet distance. Consider a man is 
walking with a dog on a leash. The man is walking on the one curve and the 
dog on the other one. Both may vary their speed, as a and j3 are arbitrary. 
Backtracking is not allowed which implies that a and f3 are non-decreasing. 
Then, the Frechet distance between the curves is the minimal length of a 
leash that is necessary. 

The Hausdorff distance between C and C rn is the smallest value d such 
that every point of C has a point of C m within distance d and every point of 
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C m has a point of C within distance d HU. It is formally defined as 
d H (C,C m ) = max{ sup inf ||a; m (£) - x(rj)\\, 

V&[^0,Vp\ ? e [°d] 

sup inf ||® m (£) -x{rj )|| }. 

£e[o,i] ■new,Vp] 

Not only (f#(C,C m ) < dp(C,C m ), but the Frechet distance between 
two curves can be arbitrarily larger than their Hausdorff distance. The 
Frechet distance is usually considered as a more reliable measure of similar¬ 
ity between curves. Figure [2] shows two curves that can be made arbitrary 
“Hausdorff-close” while being quite dissimilar. 



Figure 2: Two curves that can be made arbitrary close when e —> 0 in terms 
of their Hausdorff distance (that is exactly 3e) but not in terms of their 
Frechet distance that remains finite and equal to the diagonal of the square. 


The definition of the Hausdorff and Frechet measures, that involves in- 
fima and suprema over infinite sets of parametrizations, makes it difficult to 
devise algorithms for computing these distances between arbitrary curves. 
However, an alternative that can lead to practical algorithms is to calculate 
the Hausdorff and Frechet distances between polygonal approximations of 
the curves under consideration, as explained in the next Section. 

2.2 Distance Between Polygonal Curves 
2.2.1 Optimal Sampling of Curves 

Let us first consider the problem of approximating an arbitrary curve by a 
polygonal curve. In order to maximize the efficiency of the distance compu¬ 
tation, it is necessary to find a polygonal curve that contains as few vertices 
as possible and still approximates the original curve with sufficient accuracy. 
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Assume m + 1 points = x m (£j), i 6 [0, m], that are sampled on C m . 
This defines a polygonal curve M formed of m segments for which segment 
i goes from p* to pj + i (see Fig. [3]). Let us do the same with C and define a 
polygonal curve N composed of n + 1 points q*, i £ [0, n]. 



Figure 3: A model edge C, a mesh edge C rn and a polygonal curve M (in 
blue). 

The goal is to find sampling points on the original curve C m (resp. C ) in 
such a way that the distance between the polygonal curve M (resp. N) and 
C m (resp. C) is smaller than a threshold distance e. 

A parametrization of C m is given by ([Tj) and a parametrization of C 
is usually available from the CAD model. A possible sampling strategy for 
these curves consists in starting with the segment [*o, x p ] and refining recur¬ 
sively the discretization with points distributed uniformly in the parameter 
space, until the desired accuracy is reached. However, C is often a Bezier 
or a rational Bezier spline. As a polynomial curve, C m is also a particular 
case of a Bezier curve. For such curves, de Casteljau’s algorithm provides 
a more efficient way to refine the discretization and control the geometrical 
accuracy at the same time. 

Consider the Bernstein basis polynomials of degree p: 

B k ) tt) = (X)(i-O p - k Z k (£€[0,1]; k = 0,...,p) (2) 

where (£) = p \^ p L ^t is the binomial coefficient. Since Lagrange and Bern¬ 
stein polynomials span the same function space, we can re-write 0 as a 
Bezier curve 

*.(0 = E5? I (f)*i. €e [0,l] (3) 

i=0 

where the x^s are the control points of the Bezier curve, that form a control 
polygon (see Figure 0). The control points x^s can be computed from the 
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(v) 

node locations x^s by means of a transformation matrix 


(p) 

B^C ~ 


4 P) (£o) . 


4 p) (£0 • 

• • 4 p) (£i) 

4 p) fe) ■ 

4 p) (£p) 


A classical way to optimally sample a Bezier curve is to use de Castel- 
jau’s algorithm. A first approximation of the Bezier curve is constructec 
the single line segment between Xq and x b p (red line segment in Figure 
The distance d between this single segment and the control polygon is an 
upper bound of the distance between the curve and the segment because of 
the convex hull property. If needed, the curve is then split into two sub¬ 
curves using de Casteljau’s algorithm, each of them coinciding exactly with 
the original curve. This argument is applied recursively (see Figure [5]) to 
every sub-curve where the distance between the control polygon and the cor¬ 
responding segment is less than e. The extremities of the sub-curves at the 
finest level are thus the vertices p* (resp. q,;) of the polygonal approximation 
for C m (resp. C ). 




Figure 4: A cubic Bezier curve, its control polygon (dashed lines) and the 
coarsest approximation of the curve in red. 

A similar algorithm can be applied to rational Bezier splines. Most of 
the CAD entities can be casted as rational Bezier splines so that this optimal 
subdivision can be applied to most of the curves that are present in CAD 
models. As an example, Figure [6] presents the subdivision of a boundary 
using different values of the threshold parameter e. For the unusual cases 
involving non-standard parametrizations, the generic recursive sampling al¬ 
gorithm can be used. 
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Figure 5: A first application of de Casteljau’s algorithm. The polygonal 
approximation in red gets closer to the control polygons (and thus to the 
Bezier curves). 


2.2.2 Discrete distances between polygonal curves 

The simplest way of approximating the distances dn(C,C m ) and d-p{C,C m ) 
is to compute the discrete Hausdorff and Frechet distances 8h(M,N ) and 
5f(M,N), i.e. the Hausdorff and Frechet distances restricted to discrete 
point sets. 

Computing 8h(M, N ) consist essentially in computing both Voronoi dia¬ 
grams of the pj’s and of the q, ’s and to find Voronoi cells of pj that contains 
qj and Voronoi cells of qj that contains pj. Computing 8h{M,N) requires 
0{m + n) log(n + m) operations. In our implementation, we do not ex¬ 
plicitly construct the Voronoi diagram. An approximate nearest neighbor 
algorithm [3] is used to locate closest points, leading to the same algorithmic 
complexity. The following result holds 


d H (M,N) < 8 h (M,N) < d H (M, N) + max{D(M), D(N)}. 
where D(M) (resp. D(N)) is the maximum distance between two successive 


points in M (resp. N ). In Section 2.2.1 the accuracy e of the polygonal rep^ 


resentation is supposed to be much smaller than the actual distance between 
the CAD and the mesh in such a way that 


\d H (M,N) - d H (C,C m )\ = 0(e). 

Consequently, the discrete distance is accurate enough if and only if 

ma x{D(M),D{N)} < e, 


which implies that the approximate nearest neighbor algorithm should not 
be applied directly to the vertices of M and N , but to sets of points sampled 
at interval e on these polygonal curves instead. Note that it is still useful to 
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Figure 6: Adaptive subdivision of a 2D boundary with e = 5 x 10 2 (88 
points), e = 10 -2 (196 points), e = 10 -3 (554 points), e = 10 -4 (858 points) 
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construct the optimal polygonal approximations M and N and to sample 
these rather than obtaining directly a dense sampling of C and C m , because 
evaluating a point of a Bezier curve is computationally much more intensive 
than evaluating a point of a line segment. Unfortunately, the number of 
points to be submitted to the approximate nearest neighbor algorithm is 
clearly a problem in our context, where a relative accuracy of 10 -6 is often 
required. 

The discrete Frechet distance S F (M,N ) [9] considers only positions of 
the “leash” where its endpoints are located at vertices of the two polygonal 
curves and never in the interior of an edge. The discrete Frechet distance 
can be computed in polynomial time, i.e. in 0(mn) operations, using a 
simple dynamic programming algorithm that is described in [9]. It is very 
difficult to find a sub-quadratic algorithm that computes the discrete Frechet 
distance. We have again 

d F (M, N) < 5 f (M , N) < d F (M, N ) + max{D(M), D(N)}. 

Computing the discrete Frechet distance may be out of reach if massive 
oversampling is applied. 

2.2.3 Direct distances between polygonal curves 

Let us now consider the real Hausdorff and Frechet distances between polyg¬ 
onal curves, that is, considering the line segments per se and not only discrete 
sets of points on them. 

The computation of the direct Hausdorff distance between two polygonal 
curves is related to the Voronoi diagram of the line segments. The distance 
can only occur at points that are either endpoints of line segments or inter¬ 
section points of the Voronoi diagram of one of the sets with a segment of 
the other. This observation leads us to the following quadratic algorithm 
for computing the direct Hausdorff distance dn(M, N ). 

• Compute the bissector of all possible pairs of segments of polygonal 
curve N. Two line segment have a bisector with up to 7 arcs (lines 
and parabolas). Store all arcs in a list. 

• Compute the intersections of each arc with M. 

• Compute the distance between those intersection points and N. The 
one sided Hausdorff distance is the maximum of those distances. 

The Voronoi diagram of line segments could be theoretically computed in 
0((n + m)log(n + m)) operations [2]. Yet, it involves the computation of 
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the whole Voronoi diagram. To our knowledge, few robust implementations 
of Voronoi diagrams of line segments exist m and no extension to higher 
dimensions than two has been proposed to date. 

It is possible to compute the Frechet distance dp(M,N) between two 
polygonal curves in O (mn log(mn)) operations |2]. The algorithm is even 
more complex that for direct Hausdorff distance. 

2.3 Geometrical error based on Taylor expansions 

In the above sections, the curves C m and C have been approximated by poly¬ 
gons, so that their relative distance can be easily and efficiently computed 
without any assumption about the curves. Another approach for simplify¬ 
ing the computation of the geometrical error is to take advantage of the fact 
that the high-order nodes x % defining the mesh edge C m are located both on 
C m and on the model curve C. A Taylor expansion of the natural param¬ 
eter in the vicinity of Xj for each curve then provides an estimation of the 
geometrical error. 

Assume a curve defined by x(t), t £ [ti, ^ 2 ] - The curvilinear abscissa s(t) 
of a point x(t) of the curve is the length of the segment defined by parameter 
range [t±,t], i.e. the length of the curve from the origin x(t\) to x(t): 



We have ds = ||x t|| dt. The arc length s(t) provides the natural parametriza- 
tion x(s) of the curve: 


x(s(t)) = x(t), te[t 1 ,t 2 \, 


with || as s || = 1, where x tS is the derivative of x with respect to s. The unit 
tangent vector to the curve is computed as 

x t 

t=x >s = T—~ 

I \ X ,t 

The curvature vector k of the curve at a point x can be defined as the 
amplitude of the variations of the unit tangent t along the curve. The 
vector k = t iS is obviously orthogonal to t because t’s amplitude is equal to 
one along s. We have 



k = x, ss = 


\ x ,t\ 


3 I x ,tt || ®,t || x ,t 


x,t • x M 

II®, til 
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It is thus possible to approximate the curve with a Taylor expansion of x 
around s o = s(to) at second order as 

s 2 

x(sq + s) = x(s 0 ) + s t| ao + — k| so + 0(s 3 ). (5) 

Applying this expansion for a mesh edge C m and the corresponding model 
curve C, the geometrical error between both curves can be estimated near 
each of their common points X{ as 


8t,l ~ II h(t m — t) 

for a linear approximation and 


°T,Q ~ 


h 2 

h(tm t) H ~(k m k) 


for a quadratic approximation. In these expressions, the unit tangent vector 
t m (resp. t) and the curvature vector k m (resp. k) is computed on C m and 
(resp. C) at point Xi, and h is proportional to a “local edge length” computed 
from the Jacobian of C m . The derivatives of x m required to compute t m and 
k m can be easily obtained from Eq. ([Tj) . For the vectors t and k related to 
the model curve C, the derivatives of x are provided by the CAD model. 
The geometrical error for the whole mesh edge C rn can then be computed as 

fr=(x> 2 ) 2 . 


2.4 A simple example 

In order to illustrate the different estimates of the geometrical error de¬ 
scribed above, we consider the simple case of a particular Lame curve: 

y = \ ~ x4 ) 1 ’ x e [°> 1 1 - 

The geometrical model is a Bezier spline representing this curve with negli¬ 
gible error (see Fig. [Tj) . The spline is parametrized with a variable u E [0,1]. 

A quadratic mesh edge is built to approximate the model curve. Its end 
vertices are fixed at the extremities of the curve (i.e. u = 0 and u = 1), and 
we explore the values of 5f, 5h and 5t for different locations of the high- 
order node on the model curve. In particular, we consider 100 locations in 
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the range u G [0.2, 0.68], where the Jacobian of the edge is positive (i.e. the 
mesh is valid). 

The evolution of the discrete Hausdorff distance Sh with the high-order 
node location is shown in Fig. [8] for different values of the accuracy threshold 
e. A value of e = 10 -3 seems to yield sufficient accuracy for this curve. In 
this particular case, the discrete Frechet distance Sf is equal to the discrete 
Hausdorff distance Sh- 

Fig. m shows the quantities Sh = Sf and St normalized by their respec¬ 
tive maximum value, so that they can be compared to each other. Although 
both curves are qualitatively similar, they do not reach a minimum for the 
same high-order node location (u = 0.617 for Sh and u = 0.573 for St), 
which can be visualized in Fig. [7] Moreover, it clearly appears that the 
Taylor-based geometrical error St is a continuously differentiable function 
of the high-order node position, while the the Hausdorff and the Frechet 
distances Sh and Sf are not differentiable everywhere, in particular at their 
minimum. 

The approximate CPU time for one distance evaluation, measured on a 
recent laptop computer, is about 3.0 • 10 -2 s for Sf, and 2.4 • 10~ 3 s for Sh 
with e = 10” 3 , whereas it is only 1.1 • 10 -6 s for St- Given its continuously 
differentiable nature and its low computational cost, the Taylor-based geo¬ 
metrical error estimate St is clearly much more appropriate than the other 
distances for an optimization procedure such as the one described in Sec. [3] 

3 Mesh optimization 

In a recent paper [[233, a technique that allows to untangle high order/curvilinear 
meshes is presented. The technique makes use of unconstrained optimiza¬ 
tion where element Jacobians are constrained to he in a prescribed range 
through moving log-barriers. The untangling procedure starts from a possi¬ 
bly invalid curvilinear mesh and moves mesh vertices with the objective of 
producing elements that all have bounded Jacobians. Bounds on Jacobians 
are computed using results of papers mm- In what follows, we extend 
the optimization procedure in charge of untangling the invalid elements in 
order to take into account the geometrical error St- 

The procedure described in Ref. [23] consists in solving a sequence of 
minimization problems, where the objective function f{xi ) is composed of 
two parts £ and J~e: 

f = £ + F e 

Here Xi is the position of node i. For a node located on a boundary, it 
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Figure 7: Geometrical model representing a Lame curve (thick blue line) 
and corresponding quadratic mesh edges minimizing 5h = $F (red thin line) 
and 6t (green thin line). 



Figure 8: Case of a Lame curve: discrete Hausdorff distance 5h = Sf for 
the quadratic mesh edge as a function of the position u of the high-order 
node on the geometrical model. 
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Figure 9: Case of a Lame curve: quantities 5h = &f and St, normalized 
by their maximum value, for the quadratic mesh edge as a function of the 
position u of the high-order node on the geometrical model. 


is possible to work with the parametric coordinate (s) of the node on the 
geometrical model entity given by the CAD model. As the scale of the 
parametric coordinate can differ significantly from the scale of the phys¬ 
ical coordinates, preconditioning may then be required for the Conjugate 
Gradient to converge properly. 

The first part £ relies on the assumption that the method is provided 
with a straight-sided mesh of high quality. This mesh has potentially been 
defined to satisfy multiple criteria, such as a predetermined size field, or 
anisotropic adaptation. The conversion of such meshes to high order is 
expected to preserve as much as possible all these features. Therefore, the 
nodes shall be kept as close as possible to their initial location in the straight 
sided mesh. In this work, the definition of £ is the one of [21] i.e. 

£{Xi) = ^jY^\\ x i~ x i\\ 2 (6) 

i 

with Xi the position of the node i in the straight-sided mesh, Kg a non 
dimensional constant and L a characteristic size of the problem. 

The second part T of the functional controls the positivity of the Jaco- 
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bian. A log barrier [21] prevent Jacobians from becoming too small: 


^(*0 = £!>*(*<, <0 


e 


with l iterating on all coefficients Bf of the Bezier expansion of the Jacobian 
of e and where 




is the log barrier function defined in such a way that T blows up when 
Bf /,/q — > e, but still vanishes when Bf = Jq. In this expression, x e is the 
vector gathering the positions Xi of all nodes i belonging to element e, and 
Jq is the constant straight sided Jacobian of e. 

The value of Kg has little influence on results. The presence of £ prevents 
the problem from being under-determined, and it orients the optimization 
procedure towards a solution that tends to preserve the straight-sided mesh, 
but it is clearly dominated by J- e when invalid elements exist in the domain. 

A Conjugate Gradient algorithm is used to minimize the objective func¬ 
tion / with respect to the node positions Xi for a fixed value of the log 
barrier parameter e. A sequence of such minimization problems is solved, in 
between which the e is progressively increased, so that the Jacobian of all 
elements is forced to exceed a user-defined target value. 

In the present work, the procedure described above is followed in a first 
step. In a second step, a similar procedure is used, albeit with an objective 
function f taking into account the geometrical model: 


f = £ + B e + V e , 


The third part T> e i of the functional controls the error 5t between the 
mesh and the geometrical model. Again, a log barrier is used: 


V e ,( Xi ) = K v ^D b ( Xi ,e) 


b 

with b iterating on all boundary mesh edges and 



where <5^ is the is the geometrical error 5t for the boundary mesh edge 6, 
and Jq is a target value for 5^. The vector x b collects the positions Xi of 
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all nodes i belonging to the boundary mesh edge b. Derivatives of St with 
respect to Xi are computed using finite differences. 

In this second step, is used as a fixed log barrier (constant e) that 
is meant to prevent the Jacobian of elements to fall back below the target 
value reached in the first step. On the contrary, T> e > is a moving log barrier 
where the parameter e' is iteratively updated to drive the geometrical error 
St towards its target value. The parameter Kx> reflects the weight given to 
the geometrical error contribution with respect to the contribution of the 
Jacobians. In this work, we typically choose values around Kjj = 0.1. 


4 Examples 
4.1 NACA0012 


We consider the classical geometry of the NACA0012 airfoil with unit chord 
length. Sequences of 6 triangular meshes are generated, where the airfoil 
is discretized by 4, 6, 10, 18, 34 and 66 elements respectively, for a total 
of 244, 298, 436, 546, 832 and 1178 elements in each mesh. The first se¬ 
quence consists of linear meshes. Two others are composed of quadratic 
meshes resulting from the optimization procedure described in Section [3j 
one is optimized for element validity only, while the other also minimizes 
the geometrical error. The last two sequences are made up of cubic meshes 
generated in the same manner as the quadratic ones. 

Figure 10 shows how the geometrical error St and the model-to-mesh 
Hausdorff distance Sh evolve with the mesh size. In meshes optimized for 
validity only, both St and Sh clearly decrease with decreasing mesh size, 
but most meshes are too coarse to yield the optimal convergence rate. How¬ 
ever, minimizing the geometrical error St reduces both St and Sh of at least 
one order of magnitude in most cases. Geometrically-optimized quadratic 
meshes are even more accurate than valid cubic meshes. The geometrical 
optimization is most beneficial for coarser meshes, which lie precisely the 
range of mesh size where refining brings less geometrical accuracy. For fine 
cubic meshes, the improvement is less significant: meshes optimized for va¬ 
lidity only converge at near asymptotic rate for both while St and Sh, while 
the convergence rate with meshes optimized for both validity and geomet¬ 
rical accuracy is not improved (or even reduced for St)- Thus, optimizing 
meshes with respect to St may be particularly interesting with numerical 
schemes of very high order running on coarse meshes, where it may yield a 
suitable geometrical approximation of the model without unnecessary mesh 
refinement. 
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Figure 10: Geometrical error St (top) and model-to mesh Hausdorff distance 
5h (bottom) for the NACA0012 profile: linear meshes (p = 1) as well as 
quadratic (p = 2) and cubic (p = 3) meshes optimized either for element 
validity only (J) or for both element validity and geometrical error (J + St)- 
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Examples of meshes are shown in Figures 11 and 12 In geometrically' 


optimized meshes, the high-order nodes located on the boundary are clearly 
moved along the CAD curve to minimize 5t, while they remain in the mid¬ 
dle of corner nodes when validity only is considered. In coarse meshes, some 
elements need to be strongly deformed to satisfy both geometrical and va¬ 
lidity criteria, which may affect the simulations adversely. Indeed, highly 
distorted elements are known to harm the accuracy of finite element ap¬ 
proximations 0. They may also deteriorate the conditioning of the spatial 
discretization operator, with negative impact on time integration [23]. More¬ 
over, a correct integration of polynomial quantities on such elements may 
require costly higher-order quadrature rules. Fortunately, the effect is less 
pronounced in finer meshes. 


In order to illustrate the impact of the geometrical accuracy on simu¬ 
lations, computations solving the Euler flow around the NACA0012 airfoil 
at Mach number M = 0.5 and 3° angle-of-attack have been carried out. 
A high-order (p = 6 i.e. sixth-order polynomials) discontinuous Galerkin 
scheme was used for the spatial discretization, and steady-state solutions 
were obtained through a pseudo-tinre approach involving a backward Euler 
scheme in combination with a Newton-Krylov solver. Slip wall conditions 
are imposed on the airfoil and characteristic-based free-stream boundary 
conditions are used at the far-held boundary of the domain. 


Results for meshes in which the airfoil is discretized with 34 elements are 
shown in Figure[l3j Unsurprisingly, the numerical method does not converge 
properly with the linear mesh, and the residual cannot be decreased by more 
than two orders of magnitude. The density held is clearly different from the 
expected solution. With the quadratic mesh optimized for validity only, the 
airfoil is represented more accurately, but the corresponding solution still ex¬ 
hibits spurious oscillations and how features near corner nodes on the wall 
boundary, where the representation of the airfoil is not smooth. A drop of 
four orders of magnitude in residual is achieved in 26 pseudo-time iterations. 
With the geometrically-optimized mesh however, the computation converges 
towards the expected smooth solution in 19 pseudo-time iterations. In this 
purely inviscid test case, the increased boundary smoothness resulting from 
the minimization of 5t is instrumental in converging towards the exact so¬ 
lution without spurious entropy generation at the boundary. Moreover, the 
geometrically-optimized quadratic mesh represents the model so accurately 
that it is meaningless to use a higher-order mesh, as the Hausdorff distance 
5h ~ 4 • 10~' 5 is probably already lower than the manufacturing tolerance of 
the airfoil. 


19 









Figure 11: Coarsest meshes of the NACA0012 profile. Top: linear mesh. 
Center: quadratic meshes optimized for validity only (left) and for validity 
as well as geometrical error (right). Bottom: cubic meshes optimized for 
validity only (left) and for validity as well as geometrical error (right). 
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Figure 12: Medium-size meshes of the NACA0012 profile. Top: linear mesh. 
Center: quadratic meshes optimized for validity only (left) and for validity 
as well as geometrical error (right). Bottom: cubic meshes optimized for 
validity only (left) and for validity as well as geometrical error (right). 
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Figure 13: Density field for the NACA0012 case: flow around the airfoil (left 
column) and zoom around the leading edge (right column). Top row: results 
with a linear mesh. Middle row: results with a quadratic mesh optimized 
for validity only. Bottom row: results with a quadratic mesh optimized for 
both validity and geometrical accuracy. 
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4.2 Rattray island 


We consider now an ocean modelling application focusing on the Rattray 
island, that is located in the Great Barrier Reef near Australia. Simula¬ 
tions have been performed, in which the shallow water equations are solved 
without diffusion nor Coriolis force. The water depth at rest is uniform and 
equal to 25 nr. Slip wall conditions are imposed on the island coast and the 
lateral sides of the domain. At the upstream and downstream sides of the 
domain, uniform free-stream conditions are prescribed with a velocity cor¬ 
responding to a Froude number of Fr = 0.02, which is representative of the 
tidal stream [25j. The island, that has a length of about 1350 m, is oriented 
at 60° compared to the free stream. In this setup, there is no source of 
vorticity, and the ideal solution is a steady irrotational flow. A view of the 
domain and the solution is given in Figure [Tf] 

A quadratic mesh of 596 elements with higher density near curved bound¬ 
aries is generated. This mesh is already valid without optimization. A sec¬ 
ond mesh is obtained by minimizing St according to the procedure described 
in Section [3j Figure [15] shows a comparison of both meshes at the tips of the 
island. The mesh boundaries of both meshes look very similar. Indeed, the 
model-to-mesh Hausdorff distance for the original mesh (Sh — 1.45 m) is not 
significantly different from the one for the optimized mesh ( Sh = 0.66 m). 
This is due to the fact that the geometry is already well resolved by the 
unoptimized mesh. However, a close examination of Figure [15] reveals that 
the boundary of the unoptimized mesh is not perfectly smooth at element 
corners, while the optimized mesh looks better in this respect. 


Even though the improvement in the Taylor-based geometrical error is 
not necessarily impressive (St = 2.19 m for the unoptimized mesh against 
St = 0.37m for the optimized one), the impact on the solution is impor¬ 
tant. Simulations were performed with the same numerical method as in 
Section 4.1 Figure 16 compares the results between both meshes. With 


the unoptimized mesh, vortices are shed from the downstream tip of the 
island, preventing the flow from reaching steady-state, while a drop of 4 
orders of magnitude in residuals can be achieved with the optimized mesh, 
leading to a nearly potential solution. This test case shows, even more than 
in the NACA0012 case, that the gain in boundary smoothness brought by 
the geometrical optimization is crucial to obtain the correct solution in some 
problems. 
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Figure 14: Rattray test case: general view of the computational domain and 
quadratic mesh (bottom left), zoom on the mesh around the island (top) and 
streamlines of the ideal solution in the vicinity of the island (right). 
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Figure 15: Rattray test case: comparison between the original quadratic 
mesh (black) and the geometrically-optimized mesh (orange) at the up¬ 
stream tip (left) and the downstream tip (right) of the island. 


4.3 High-lift airfoil 


In this section, we apply the methods described in Section [3] to an acoustic 
application involving a high-lift airfoil. The geometry is a 3-element airfoil 
based on the RA16SC1 profile, with the slat and flap deflected by 30° and 
20° respectively. The chord of the main element is 480 mm, and the compu¬ 
tational domain is a disc of radius 1 m centered on a point P located close 
to the trailing edge. The acoustic excitation consists of a monopole source 
placed at point P, with an amplitude of 1 Pa and frequency of 7816 Hz. The 


computational domain is shown in Figure 17 


Simulations are performed with a discontinuous Galerkin nodal scheme 
for the space discretization, and the standard fourth-order four-stage Runge- 
Kutta integrator in time. Slip wall conditions are imposed on the airfoil, 
while characteristic-based non-reflecting boundary conditions are prescribed 
in the far field. The monopole is modeled by a Gaussian pressure pertur¬ 
bation of half-width 3 mm. Simulation are run until a periodic regime is 
reached. A reference solution is obtained through a computation on a fine 
grid. 

Two sets of triangular meshes are generated, one composed of medium- 
size meshes (1617 elements) and one composed of coarse meshes (296 el¬ 
ements). Each set consists of a linear mesh, the corresponding quadratic 
mesh optimized for element validity only and the corresponding quadratic 
mesh optimized for both validity and geometrical accuracy. Details of the 
meshes around the slat and the leading edge of the main component, where 
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Figure 16: Rattray test case: contours of the sea surface elevation (upper 
row) and streamlines (bottom row) at t = 2100 s. Results obtained with the 
original quadratic mesh (left column) and with the geometrically-optimized 
mesh (right column). 
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Figure 17: Acoustics test case: computational domain and medium-size 
quadratic mesh (left). 
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the boundaries influence most the acoustic field, are plotted in Figure 18 


It is obvious that the coarse quadratic mesh optimized for validity only 
represents the model very poorly, while the minimization of the geometrical 
error yields a fairly accurate and smooth approximation of the airfoil. In 
the region of the computational domain shown in Figure [T8l the geometrical 
optimization decreases the discrete model-to-mesh Hausdorff distance Sh by 
a factor 7 approximately (Sh = 9.3 mm to Sh = 1.3 mm), and the geomet¬ 
rical error St drops by a factor 8 (St = 15.5mm to St = 1.8mm). Above 
all, a close examination of mesh optimized for validity only at the trailing 
edge of the slat shows that the mesh edge on the lower side of the profile 
crosses the edge representing the upper side: even though the mesh is valid 
in the finite element sense, it is physically incorrect. On the contrary, a 
simulation with the coarse quadratic mesh optimized for geometrical accu¬ 
racy can give an acceptable solution, provided that the spatial discretization 
is of sufficiently high order. In the present case, simulations with a 10 th - 
order discontinuous Galerkin scheme were slightly more dissipative than the 
reference computation on a fine mesh (see Figure [l9|). 

The difference between both quadratic medium-size meshes is less spec¬ 
tacular, as seen in Figure [18} Sh decreases by a factor 4 (2.5 mm to 0.6 mm) 
and St by a factor of 5 (3.8 mm to 0.8 mm) with geometrical optimization. 
Simulations are run with these meshes, and the RMS acoustic pressure prms 
is measured over the last 6 oscillation periods along a circle of radius 750 mm 
centered at point P. The results are expressed in terms of Sound Pressure 
Level as SPL = 20 log(pRMs/.Pref), where p re f = 2 • 10~ 5 Pa. Figure [20] shows 
that the effect of the geometrical optimization impacts significantly the ac¬ 
curacy of the sound directivity. 


5 Extension to three-dimensional meshes 

In the same manner as described in Section [2] for curves, it is possible to 
define a geometrical error for a surface mesh element S m approximating a 
model surface S. At each of interpolation point Xi where both surfaces 
coincide, a first-order estimation of the geometrical error is: 

S l T = \\h(n m - n)|| . 

where n m represents the unit normal to S m and n represents the unit normal 
to S. Here, h is proportional to the square root of a “local surface element 
area”, determined from the Jacobian of S m . The geometrical error St for 
the surface mesh element can then be computed from all the S l T in the 
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Figure 18: Acoustic test case: details of the mesh around the slat and the 
leading edge of the main component for the coarse meshes (left column) and 
the medium-size meshes (right column). Linear meshes (top row), quadratic 
meshes optimized for geometry only (middle row) and quadratic meshes 
optimized for both validity and geometrical accuracy (bottom row). 
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Figure 19: Acoustics test case: acoustic pressure field at final time for the 
reference solution (left) and the solution obtained on the coarse quadratic 
mesh optimized for geometrical accuracy (right). 


element. It is then possible to use the optimization process presented in 
Section [3] in order to obtain a better representation of the model surfaces by 
the boundary of a 3D volume mesh. 


In order to illustrate the potential of this approach for 3D meshes, we 
apply the method to the case of a wing made from an extruded NACA0012 
profile of unit chord length. Fig.|21|shows a coarse mesh of 847 tetrahedra in 
two versions, namely one optimized for validity only and one optimized for 
both validity and geometrical accuracy. As in the 2D case, the smoothness 
of the mesh boundary at the leading edge is significantly improved by the 
minimization of St, and the approximation of the airfoil seems to be more 
accurate. The geometrical error St is indeed decreased by more than an 
order of magnitude (0.16 to 0.012). 


Another example is the case of an ONERA M6 wing of chord length 

A coarse quadratic volume 


of 810 at wing root, illustrated in Figure 22 


mesh of 11851 terahedra is generated, then optimized for validity only on 
one hand, and for both validity and geometrical accuracy on the other hand. 
The geometrical error St is reduced from 40 to 15 by the geometric optimiza¬ 
tion. The representation of the leading edge is clearly improved, particularly 
where it merges with the tip surface of the wing. 
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Figure 20: Acoustics test case: sound directivity expressed in terms of Sound 
Pressure Level (in dB) for the reference solution, as well as the solutions 
obtained with quadratic medium-size meshes optimized for validity only (J) 
and for both validity and geometrical accuracy (J + 5t)- 
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Figure 21: Mesh of the 3D NACA0012 geometry optimized for validity only 
(left) and for both validity and geometrical accuracy (right). 

6 Conclusions 

In this paper, we have presented methods to evaluate and improve the geo¬ 
metrical representation of CAD models in high-order meshes. 

The interest of formal distances in the plane for this purpose has been 
assessed. The Frechet and the Hausdorff distances between two curves cor¬ 
responding respectively to the model and to the mesh boundary have been 
examined. A discrete version of these quantities, in particular the Hausdorff 
distance, can be computed fast enough to assess the quality of the geo¬ 
metrical model approximation in practical 2D meshes. However, it is still 
computationally too costly to be employed in mesh optimization algorithms. 

To this end, a fast estimate of the geometrical error between the mesh 
boundary and the CAD model is presented, which is based on a Taylor 
development of each curve. It is then introduced in a pre-existing optimiza¬ 
tion framework that guarantees the mesh validity. Several examples show 
that minimizing this quantity significantly improves the representation of 
the model: depending on the case, the model-to-mesh Hausdorff distance 
is decreased by a factor 2 to 35, the gain being larger for coarse meshes. 
An important aspect of the method lies in the beneficial impact of the ge¬ 
ometrical optimization on the mesh boundary smoothness. As evidenced 
by several test cases, this effect is often instrumental in obtaining accurate 
solutions from high-order simulations. The approach is easily extended to 
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Figure 22: ONERA M6 case: General view of the volume mesh and the wing 
(top), details of the wing surface mesh optimized for validity only (center 
left and bottom left) and for both validity and geometrical accuracy (center 
right and bottom right). 
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3D meshes, as illustrated by two examples. 

The method presented in this paper reduces the need to refine a high- 
order mesh only for the purpose of representing the geometrical model cor¬ 
rectly. Therefore, it makes it easier to enjoy the computational efficiency 
of very high order numerical schemes in practical simulations. However, 
the constraints of element validity and geometrical accuracy imposed on the 
mesh may lead to highly curved elements inside the computational domain. 
The impact of the element distortion on the accuracy, computational cost 
and robustness of the simulation remains to be assessed and, if possible, 
controlled. This topic is the subject of ongoing work. 
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